Reparameterization of GFN1-xTB for atmospheric molecular clusters: applications to multi-acid–multi-base systems

Atmospheric molecular clusters, the onset of secondary aerosol formation, are a major part of the current uncertainty in modern climate models. Quantum chemical (QC) methods are usually employed in a funneling approach to identify the lowest free energy cluster structures. However, the funneling approach highly depends on the accuracy of low-cost methods to ensure that important low-lying minima are not missed. Here we present a reparameterized GFN1-xTB model based on the clusteromics I–V datasets for studying atmospheric molecular clusters (AMC), denoted AMC-xTB. The AMC-xTB model reduces the mean of electronic binding energy errors from 7–11.8 kcal mol−1 to roughly 0 kcal mol−1 and the root mean square deviation from 7.6–12.3 kcal mol−1 to 0.81–1.45 kcal mol−1. In addition, the minimum structures obtained with AMC-xTB are closer to the ωB97X-D/6-31++G(d,p) level of theory compared to GFN1-xTB. We employ the new parameterization in two new configurational sampling workflows that include an additional meta-dynamics sampling step using CREST with the AMC-xTB model. The first workflow, denoted the “independent workflow”, is a commonly used funneling approach with an additional CREST step, and the second, the “improvement workflow”, is where the best configuration currently known in the literature is improved with a CREST + AMC-xTB step. Testing the new workflow we find configurations lower in free energy for all the literature clusters with the largest improvement being up to 21 kcal mol−1. Lastly, by employing the improvement workflow we massively screened 288 new multi-acid–multi-base clusters containing up to 8 different species. For these new multi-acid–multi-base cluster systems we observe that the improvement workflow finds configurations lower in free energy for 245 out of 288 (85.1%) cluster structures. Most of the improvements are within 2 kcal mol−1, but we see improvements up to 8.3 kcal mol−1. Hence, we can recommend this new workflow based on the AMC-xTB model for future studies on atmospheric molecular clusters.


Introduction
Molecular clusters, formed through the aggregation of various atmospheric species, play a central role in aerosol particle formation. 1Aerosols are liquid or solid ne particles suspended in air that can act as cloud condensation nuclei (CCN) if they reach sizes at or above 50-100 nm. 2 Roughly 50% of CCN are believed to be initially formed as clusters. 3CCN acts as nucleation cores for water uptake and then further growth into clouds meaning there is a direct correlation between aerosols and cloud number/properties and hence the climate.The biggest uncertainty in modern climate forcing predictions is due to the uncertainties from aerosol-cloud interactions. 1ulfuric acid has been shown to be the main driver of cluster formation.5][6][7] It is extremely difficult to experimentally measure the composition and formation mechanism of the initial clusters due to their small size and neutral charge.][10] This leaves theoretical studies as the only way to elucidate the thermodynamics, kinetics, and molecular interactions governing cluster formation and its evolution.The main challenge for studying atmospheric molecular clusters is their complex congurational spaces, which require advanced congurational sampling techniques and computationally demanding quantum chemistry methods to evaluate the cluster properties accurately. 5Furthermore, atmospheric clustering is believed to be a multi-species process, 11 adding another dimension of chemical complexity.
Thoroughly exploring the congurational space of atmospheric molecular clusters using, for instance, metadynamics simulations 12,13 or genetic algorithms [14][15][16] at a high level of theory is extremely computationally demanding.Hence, usually, a funneling approach 5,7,17 is applied, where the congurational space is initially explored at a low level of theory such as force-eld or semiempirical methods, and only a subset of low energy structures is selected, reoptimized, and reexamined at a higher level of theory.This process is repeated with an increasing level of theory until only a few structures remain for evaluation at the desired high level.Schematically, the process can be given as: (1) Generate initial cluster congurations: ABCluster/OGOLEM/Basin hopping or similar.(2) Semi-empirical calculations: Optimization at the PM6/PM7/GFN1-xTB/GFN2-xTB or similar level.
(3) DFT calculations: DFT optimization and vibrational frequency calculations.(4) Single point energy renement: Single-point energy calculation at coupled cluster level on the DFT optimized geometry.
Between each step in the funneling approach, ltering can be applied to reduce the number of structures that need to be handled.This can either be based on an energy threshold or a set number of cluster structures.Eventually, we end up with a handful of structures at the highest obtainable level.
The rst step in the funneling procedure is the generation of a large number of conguration candidates.The key idea is sampling a large part of the potential energy surfaces at a low level of theory to get estimates for the global free energy minimum.9][30][31][32] The major issue at this step is that most force-eld methods are unable to describe bondbreaking, such as proton transfer reactions, which are important for atmospherically relevant molecular clusters, requiring the sampling to include monomers where the hydrogens have been transferred to get adequate sampling.Furthermore, the accuracy of force-eld methods is also insufficient to determine a subsample of the conformer candidates and all the candidates have to be taken to the higher level of theory.
The next step is semi-empirical calculations as these are a better description of the chemistry and ltering can be applied.6][7] GFN stands for geometries, frequencies and non-covalent interactions, which are the main target properties for the method.PM stands for parameterization method indicating the model version.Of these methods, GFN1-xTB has shown to have the highest correlation with electronic binding energies at a higher level of theories 37,38 and have been shown to have a higher correlation with DFT trajectories for molecular dynamics than GFN2-xTB. 39The reason GFN2-xTB performs worse than GFN1-xTB for atmospheric molecular clusters (oen involving sulfuric acid) is that there is a decrease in the number of d-functions for sulfur in the basis set for the newer GFN2-xTB model.
The third step is the subsequent optimization and vibrational frequency calculation of the structures with DFT.This is the main bottleneck in the sampling methodology as limited computational resources only allow a xed number of DFT structures to be optimized.Therefore some form of ltering is required, oen based on structural properties or electronic energies from the semiempirical calculations.To circumvent the inaccuracies of semiempirical methods an intermediate step can be included, involving single-point energy calculations at the DFT level on as many structures as possible.Another option for an intermediate step is the utilization of machine learning (ML) methods.One can calculate a subset of the structures at a desired DFT level and train an ML model to predict the energies of the remaining structures. 39,40However, to mimic accurate DFT energies, kernel-based ML methods become computationally demanding [40][41][42] and neural-networks will require an extensive set of training data and hyperparameter optimization. 43Moreover, ML methods oen fail when predicting on structures different from the training set.
Overall, the funneling approach is never more efficient than its weakest link given by the semiempirical step in 2, in which accuracy determines the number of structures that have to be optimized/have single points calculated at the DFT level.In this paper, we focus on reparameterizing the GFN1-xTB method based on DFT energies of atmospherically relevant molecular clusters yielding a GFN1-xTB model reparametrized based on uB97X-D/6-31++G(d,p) for 'atmospheric molecular clusters' denoted AMC-xTB.This new parameterization is used to sample 288 large multi-acid-multi-base clusters containing AM equivalent to the clusters studied by Knattrup et al. 44 2 Methodology

Computational details
Single-point energies, gradients, and geometries for the reparameterization, congurational sampling, and comparisons were calculated using the xtb 6.4.0 program using the GFN1-xTB 33 and AMC-xTB parameterizations.A modied version of ArbAlign 45 available in the JKCS program 46 was used to calculate the root-mean-square differences (RMSD) between molecular structures.Gaussian 16, version B.01 47 was used for the DFT calculations.CREST 2.12 12,13 with an energy window of 15 kcal mol −1 and in noncovalent interaction mode and ABCluster 2.0 15,16 with a population of SN = 3000, number of generations of g max = 200, and gen.survival of g limit = 4 were used for additional congurational sampling.

Cluster data sets
For reparameterization of GFN1-xTB, we used the clusteromics I-V data sets [48][49][50][51][52] containing (acid) 0-2 (base) 0-2 clusters of the following atmospherically relevant species: sulfuric acid (SA), methanesulfonic acid (MSA), nitric acid (NA), formic acid (FA), ammonia (AM), methylamine (MA), dimethylamine (DMA), trimethylamine (TMA) and ethylenediamine (EDA).All structures are optimized at the uB97X-D/6-31++G(d,p) level of theory, as benchmark studies 37,53 show this to be a good compromise between accuracy and speed, and we used up to the three lowest electronic energy congurations found per each cluster as the optimization set for GFN1-xTB reparametrization.This leads to an optimization set comprising of a total of 1073 clusters.The GFN1-xTB reparameterization based on this optimization set will be denoted as the AMC-xTB model.
All new data calculated is freely available in the Atmospheric Cluster DataBase 54 along with the new AMC-xTB parameter le (see Section S2).†

Optimization strategy
The GFN1-xTB model contains 15 global parameters, 2 element-pair-specic parameters, and 32 element-specic parameters of relevance to the species present in the optimization sets (H, C, N, O and S atoms).Initially, the Hessian was generated to probe the sensitivity of the parameters, however, we found it computationally feasible to employ a similar optimization strategy to the original GFN1-xTB paper, 33 where we optimize all relevant parameters simultaneously.We utilize a modied version of an in-house pseudo-Newton-Raphson optimizer by Jensen et al. 55 for the optimization of a target function (T) containing a linear combination of the difference in electronic binding energies (DE b ) in kcal mol −1 and gradient norms (g norm ) in hartree bohr −1 radius at the current GFN1-xTB parameterization and the target uB97X-D/6-31++G(d,p) level of theory: here, N conf is the total number of structures in the optimization set, N atoms i is the number of atoms in the i-th structure.We normalized by the number of atoms in each structure to prevent "overtting" to the larger clusters.
We chose the electronic binding energies as the target properties to get a better tool for ltering based on energies in congurational sampling procedures.The gradients were included directly in the target function.We use equilibrium structures at the given level of theory, which are supposed to have near-zero gradients.However, the upper limit is set by the accuracy threshold within the xTB program during the optimization, which makes the gradients non-zero in the calculations.
Including only the electronic bindings energies in the target function yields a much better t for the energies but causes the gradients to "explode", effectively rendering the optimization functionality of AMC-xTB useless for our target species.Giving higher weight to the gradient norms in the target function makes the optimized structures more similar to the target uB97X-D/6-31++G(d,p) level of theory, however, we found that it causes problems in the congurational sampling procedure where the decreased accuracy in determining the binding energies causes our congurational sampling to yield high-energy conformers at the DFT level.Overall, we found that including the gradient norms and difference in electronic bindings energies in a 1 : 1 ratio as the best compromise between the two properties.

Updated congurational sampling workows
The strength of the new AMC-xTB model is that it can be used directly in congurational sampling programs such as ABCluster and CREST.Here, we will test two different new workows for applying the reparameterized models in cluster congurational sampling.
2.4.1 Original workow.The workow usually employed in studying atmospheric molecular clusters can be summarized as: here the number of congurations N that have to be optimized at the DFT level is a severe bottleneck in the number of new cluster systems that can be studied.Usually, 50-100 congurations are optimized at the DFT level.2.4.2Independent workow.The independent workow refers to congurational sampling from scratch using the wellestablished funneling approach using AMC-xTB instead of GFN1-xTB. 5,17As the aim for the approach is to be generally applicable, we also included an additional CREST step, as it should be better at handling exible organic compounds: here, ABCluster, a genetic algorithm for sampling clusters, is used for the initial sampling of all possible neutral/ionic combinations of monomers that yield overall neutral clusters.The xTB 6.4.0 program was then used to optimize all the congurations at the AMC-xTB level.The cluster lowest in electronic energy was then taken as the input structure for CREST in non-covalent interaction mode, again using our new AMC-xTB model.The initial ABCluster sampling is needed because we found CREST to be quite sensitive to the input structure and, therefore, needs a good guess for a starting structure.The 50 structures lowest in electronic energy are then optimized at the DFT level.2.4.3Improvement workow.The improvement workow refers to using the best structure currently known at the corresponding level of theory as the input structure for CREST using the AMC-xTB model.
Best structure / CREST(AMC-xTB) / DFT opt 50 lowest from here on, the workow is the same as the independent workow, where the 50 structures lowest in electronic energy are optimized at the corresponding DFT level.

Extension of the multi-acid-multi-base dataset
To gain a more complete test set we extended the multi-acidmulti-base clusters systems by Knattrup et al. 44 using the same workow for a total of 288 new AM-containing clusters.With the acids being SA, MSA, FA and NA and the bases being AM, MA, DMA and TMA.This is the rst sampling of multi-component clusters containing up to 8 different species yielding a data set where synergistic effects in cluster formation between different species of bases 41 and acids 44 and mixes thereof can be studied.Such clusters could be relevant for modeling polluted coastal environments.Fig. 1 presents the molecular structure of a newly sampled 8-component cluster.It is seen that all the acids have transferred a proton to all the bases.

Assessment of the AMC-xTB binding energies
We reparameterized GFN1-xTB to obtain a new tight-binding semiempirical reparameterization denoted as AMC-xTB.Fig. 2 shows the error in electronic binding energies before (GFN1-xTB) and aer (AMC-xTB) reparameterization.The models have been tested on the entire clusteromics I-V [48][49][50][51][52] data sets (56 436 data points), the sulfuric acid-multi-base (SA) 1-4 (AM/MA/ DMA/TMA/EDA) 1-4 cluster data set (684 data points) by Kubečka et al. 41 and the multi-acid-muti-base (SA/FA/MSA/NA) 1-4 (MA/ DMA/TMA) 1-4 by Knattrup et al. 44 including the new AMcontaining clusters (1629 data points).All the tested structures are equilibrium structures at the uB97X-D/6-31++G(d,p) level of theory.Although the Gaussian version and integration grid used for optimization differ for some structures, it was found to have a negligible effect on this comparison as we are studying the binding energies and not the absolute energies.
For all the data sets shown in Fig. 2 the reparameterization results in near-zero means of the energy errors.This is a reduction from error means of 3.7-11.8kcal mol −1 for GFN1-xTB.In addition, the AMC-xTB model achieves a more narrow error distribution with the root mean square deviations decreasing from 7.6-12.3kcal mol −1 to 0.81-1.45kcal mol −1 , implying that there will be fewer outliers.The error span on the larger clusters for the Knattrup et al. 44 and Kubečka et al. 41 sets are similar to the error span for the smaller clusters in the optimization set.This shows that reparameterizing on smaller clusters is adequate for calculations on larger-sized clusters as the model gets some of the underlying chemistry correct and scales effectively with system size.The new AMC-xTB model reduces the number of structures needed to pass from the semiempirical step to the DFT step in congurational sampling.For atmospheric molecular clusters, this implies that the AMC-xTB model is unequivocally better to apply in the congurational sampling funneling workow compared to GFN1-xTB.

Assessment of the AMC-xTB geometries and gradients
The gradient norms were also a part of the optimization scheme (eqn (1)).Fig. 3 shows the gradient norms given by the xtb program for the two parameterizations.The structures are equilibrium structures at the uB97X-D/6-31++G(d,p) level of theory, so the ideal gradient norms should be below the default gradient convergence thresholds of 10 −3 E h /a.None of the methods manages to be below this threshold, but the AMC-xTB model is close.This does not directly mean the model is closer to the correct structure, as the new parameters might just have attened the potential energy surface at this point without moving closer to the minimum.However, including the gradients in the target function avoids numerical instability, and yields reasonable optimized structures.To test if the structures are closer to a minimum at the DFT level, the initial DFT structures of all three datasets were optimized using the different parameterizations, and the RMSD was computed between the initial DFT structure and the GFN1-xTB/AMC-xTB optimized structures (see Table 1).
We nd that the AMC-xTB model reduces the mean RMSD of the full clusteromics set from 0.484 Å to 0.355 Å and a similar reduction is seen for the Knattrup et al. 44 and Kubečka et al. 41 Fig. 1 The (SA) 1 (MSA) 1 (NA) 1 (FA) 1 (AM) 1 (MA) 1 (DMA) 1 (TMA) 1 cluster structure lowest in Gibbs free energy at the uB97X-D/6-31G++(d,p) level of theory with the quasi-harmonic approximation (cutoff of 100 cm −1 ) for the initial sampling.Yellow = sulfur, red = oxygen, cyan = nitrogen, brown = carbon and white = hydrogen.sets with RMSDs being reduced from 0.378 Å to 0.235 Å and from 0.376 Å to 0.189 Å, respectively.This, coupled with the smaller gradients, suggests that the reparameterized model is closer to a minimum at the DFT level.This implies that the preoptimization step in a funneling approach with the AMC-xTB model compared to GFN1-xTB yields structures closer to the DFT structure and will likely reduce the subsequent optimization time at the DFT level.

Test of new congurational sampling workows
To further test how the new AMC-xTB model fares in cluster congurational sampling, we tested the independent and improvement workow for several previously studied (acid) 4 (base) 4 cluster systems.Hence, the workow is tested on clusters up to twice the size of those used in the reparameterization.
In the case of the (SA) 4 (AM) 4 and (SA) 1 (MSA) 1 (NA) 1 (FA) 1 (-AM) 1 (MA) 1 (DMA) 1 (TMA) 1  clusters, the independent/ improvement workows perform similar and yield similar free energy improvements.However, for the (SA) 4 (AM) 1 (MA) 1 (-DMA) 1 (TMA) 1 clusters, the independent workow works slightly better, nding a cluster 0.85 kcal mol −1 lower in free energy compared to the improvement workow.This illustrates that the sampling is very sensitive to the conguration used as input for CREST, although it might also be due to the randomness of the dynamic processes in CREST.The reason for the difference might be that the original work's congurational sampling was worse than the independent workow, yielding a worse starting structure for the CREST sampling within the improvement workow.We see a massive improvement in the congurational sampling of the (SA) 4 (EDA) 4 clusters by −18/−21 kcal mol −1 .This is caused by the exibility of the EDA molecule, as it is the only monomer that contains a C-C bond it can rotate around, making it difficult to sample the full congurational space using only ABCluster with rigid molecules.This improvement should primarily be attributed to the inclusion of metadynamics sampling in CREST and not purely the parameterization of AMC-xTB as it allows the EDA to rotate around its bonds and nd a structure with more/better paired intermolecular interactions as seen in Fig. 5.It should also be noted that the main Fig. 3 The gradient norms for the GFN1-xTB and AMC-xTB methods.a is the Bohr radius.The clusteromics [48][49][50][51][52] set is (SA/FA/MSA/NA) 0- Table 1 Comparison of the mean, median, and standard deviation (std) of the root-mean-squared differences (RMSD) between the initial DFT structure and the optimized structure at the given parameterization.The clusteromics I-V [48][49][50][51][52] sets includes the (SA/FA/MSA/NA) 0- 2 (AM/MA/DMA/TMA/EDA) 0-2 clusters, the Kube čka et al. 41 set comprise the sulfuric acid-multi-base (SA) 1-4 (AM/MA/DMA/TMA/ EDA) 1-4 clusters and Knattrup et al. 44 set has the multi-acid-muti-base (SA/FA/MSA/NA) 1-4 (MA/DMA/TMA)  improvements are the electronic binding energy and the thermal contribution varies very little between the clusters.However, this shows the strength of the presented workows as they can be used for clusters containing more exible organic molecules.

Massive improvement test
Based on the previous sections, it is clear that the improvement workow could locate more stable clusters.As the potential energy surface of multi-acid-multi-base clusters becomes very complicated, here we test this new approach for such systems.These 288 new AM-containing clusters were used as a massive test set for the improvement workow using the newly parameterized AMC-xTB model.The improvement workow manages to nd congurations lower in free energy at the uB97X-D/6-31++G(d,p) for 245 out of the 288 clusters (85.1%) as can be seen in Fig. 6.
Comparing the conformer index at the AMC-xTB level of theory with the nal uB97X-D/6-31++G(d,p) level of theory with the quasiharmonic approximation (cutoff of 100 cm −1 ) the Gibbs free energy minimum at the DFT level is also the electronic energy minimum at the AMC-xTB level of theory for 66 of the clusters (see Fig. S1 †).If 10 conformers are included from the AMC-xTB level of theory, the free energy minimum energy is captured for 155 out of the 288 clusters with improvements found for 209 (see Fig. S2 †).Furthermore, the maximum error is 2 kcal mol −1 with a mean of 0.12 kcal mol −1 when reducing from 50 to 10 conformers.
This highlights the need for including dynamics-based sampling procedures for atmospheric clusters even though the system might seem fairly rigid.It can also be envisioned that the improvement workow will be quite important when studying much larger (SA) 1-20 (base) 1-20 clusters as recently done by Engsvang et al. 42,62 and Wu et al. 38 For large clusters, the global minimum is tricky to locate, and adding dynamics-based congurational sampling might aid in the process.

Conclusions
We have reparameterized the GFN1-xTB model to yield better binding electronic energies and gradient norms for atmospherically relevant clusters composed of the following species: sulfuric acid (SA), methanesulfonic acid (MSA), nitric acid (NA), formic acid (FA), ammonia (AM), methylamine (MA), dimethylamine (DMA), trimethylamine (TMA), and ethylenediamine  (EDA).The reparameterization, denoted AMC-xTB, for use in the xtb/CREST program, is based on the uB97X-D/6-31++G(d,p) level of theory.The model shows a substantial decrease in the error of the binding electronic energies compared to the original GFN1-xTB parameterization and the gradient norms of the equilibrium structures are closer uB97X-D/6-31++G(d,p) level of theory compared to GFN1-xTB.The reparameterization strategy is general and can be used to reparameterize other methods such as GFN2-xTB.
We tested two new congurational sampling procedures with the new parameterizations being employed in the xTB and CREST programs.The rst workow, denoted as "improvement workow," is based on improving the best structure currently known in the literature with CREST and then doing the DFT calculations.The second workow, denoted the "independent workow," starts by congurational sampling using ABCluster, followed by xtb, CREST, and then DFT.Using the two workows, we nd cluster structures lower in free energy for the following (SA) 4 (EDA) 4 , (SA) 4 (AM) 4 ,(SA) 4 (AM) 1 (MA) 1 (DMA) 1 (TMA) 1 , (SA) 1 (-MSA) 1 (NA) 1 (FA) 1 (AM) 4 and (SA) 1 (MSA) 1 (NA) 1 (FA) 1 (AM) 1 (MA) 1 (-DMA) 1 (TMA) 1 systems in all cases compared to the best-known value in the literature.
Testing the improvement workow on 288 large multi-acidmulti-base cluster systems, the workow nds improvements for 85.1% of the clusters, showing the need for dynamics-based sampling.
The parameterization strategy given here is not specic to either GFN1-xTB or atmospheric clusters and could be used in general.For instance, one could imagine increasing the number of d-functions in the basis set for sulfur atoms in GFN2-xTB and then reparameterizing the new GFN2-xTB model or doing a reparameterization for much larger clusters.

Fig. 4
Fig.4Comparison of the lowest free energy conformer found by the independent and improvement configurational workflows compared to the configurations found by Kube čka et al.41 (a-c) and two new multi-component AM clusters sampled in the same way as Knattrup et al.44 (d and e).Gibbs free energies are calculated at the uB97X-D/6-31++G(d,p) level of theory with the quasi-harmonic approximation (cutoff of 100 cm −1 ) and vib.frequencies scaled by 0.996 in accordance with Kube čka et al.41

Fig. 6
Fig.6Comparison of the lowest free energy conformer found by the improvement configurational workflow compared to the original workflow.Gibbs free energies are calculated at the uB97X-D/6-31++G(d,p) level of theory with the quasi-harmonic approximation (cutoff of 100 cm −1 ).
1-4 clusters including the new AMcontaining clusters sampled in this work.The lowest errors are shown in bold